function [gradX,gradY]=grad_geo(lon,lat,data)

% Function to compute meridional and zonal gradients
% in the case of an irregular grid
%
% Input: 
%       - lon : longitude in degrees
%       - lat : latitude in degrees
%       - data : data(latitude dimension,longitude dimension)
%
% G. Charria 12/2007

gradX(1:size(data,1),1:size(data,2))=NaN;
gradY(1:size(data,1),1:size(data,2))=NaN;

%1D case
if ((size(lon,1)==1) | (size(lon,2)==1))
    for ilon=2:length(lon)-1
        for ilat=2:length(lat)-1
            gradX(ilon,ilat)=(data(ilon+1,ilat)-data(ilon-1,ilat))./((lon(ilon+1)-lon(ilon-1)).*2.*pi.*6378e3/360.*cos(lat(ilat).*pi./180));
            gradY(ilon,ilat)=(data(ilon,ilat+1)-data(ilon,ilat-1))./((lat(ilat+1)-lat(ilat-1)).*2.*pi.*6378e3./360);
        end
    end
else % 2D case
    for ilat=2:size(lat,1)-1
        for ilon=2:size(lon,2)-1
            gradX(ilat,ilon)=(data(ilat,ilon+1)-data(ilat,ilon-1))./((lon(ilat,ilon+1)-lon(ilat,ilon-1)).*2.*pi.*6378e3./360.*cos(lat(ilat,ilon).*pi./180));
            gradY(ilat,ilon)=(data(ilat+1,ilon)-data(ilat-1,ilon))./((lat(ilat+1,ilon)-lat(ilat-1,ilon)).*2.*pi.*6378e3./360);
        end
    end
end
    
    
    

